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SUMMARY 

The three-dimensional marginal separation of a boundary layer along a line of 
symmetry is considered. The key equation governing the displacement function is derived, 
and found to be a non-linear integral equation in two space variables. This is solved 
iteratively using a pseudo - spectral approach, based partly in double Fourier space, and 
partly in physical space. Qualitatively the results are similar to previously reported 
two-dimensional results (which are also computed to test the accuracy of the numerical 
scheme); however quantitatively the three-dimensional results are much different. 


♦Work funded under Space Act Agreement C99066G; presently at Department of Mathematics, University 
of Manchester, Manchester, England. 
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1. INTRODUCTION 


High Reynolds number flows may be broadly categorised into two classes, namely 
attached and detached. The former class generally involves the confinement of a viscous 
(i.e. boundary) layer to the immediate vicinity of the flow boundary, the flow elsewhere 
being generally inviscid. A classical example of this is the Blasius flow over a 
semi— infinite flat plate. The detached class of high Reynolds number flow involves 
viscous layers occuring in the main body of the fluid (as well as on the flow boundaries). 
An important example of a flow of this class is that of the steady flow past a circular 
cylinder. This problem has been investigated numerically (for finite Reynolds numbers) by 
Fomberg (1980,1985) and models for asymptotically large Reynolds number situations been 
described by Smith (1979,1985) and by Peregrine (1985), although the precise details of 
the ultimate high Reynolds number flow in this, (and similar important cases) are still 
open to speculation. Generally, attached flows may be thought of as being associated with 
streamlined bodies, whilst detached flows are more likely to be associated with bluff 
bodies. 

One reason for the great theoretical difficulty in describing detached flows is that 
generally there appears to be no continuous transition from the fully attached to the 
grossly detached state. A classical example of this is the work of Goldstein (1948), who 
showed convincingly, for the first time bow a boundary layer with an adverse pressure 

gradient would terminate in a singularity right at the the point of separation. Some years 

later Stewartson (1970) showed that this singularity could not be removed or alleviated 
with the inclusion of a triple deck. 

The problem studied by Stewartson et al (1982) (hereafter refered to as I) was that 
of two-dimensional marginal separation. Here the question was the behaviour of a 

boundary layer which is on the verge of separating at a point, but then recovers. This 

work showed how an interaction region was formed, centred at the point of zero wall 
shear, which served to smooth - out the discontinuous stream wise gradient of the shear 
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stress at this point. Further, it was shown in I (by the solution of a non-linear integral 
equation) that by means of varying an arbitary parameter 

( T in their notation) that : (i) separated flow was possible within the interaction zone; 
(ii) for certain values of P the solution exhibited non— uniqueness; (iii) there exists a 
critical value of P above which no solution exists. (This latter conclusion is a further 
example of a lack of uniform transition from a fully attached to a grossly detached state.) 
A thorough investigation of the aformentioned integral equation was mounted by Brown 
and Stewartson (1983), and this revealed that up to four solutions were possible for certain 
values of the parameter P. 

One of the motivations for this work was some numerical investigations of laminar 
boundary layers on smooth bodies of revolution at incidence, carried out by Cebeci et al 
(1980). In this study, the classical boundary layer equations along the line of symmetry 
were integrated numerically, and it was found that at one particular angle of incidence 
(40° for the case cited), the skin friction dipped (linearly) to zero at one location, and 
then rose quickly again. At a very slightly increased angle of incidence (41 e ) the 

separation process turns out to be catastrophic (although this is not unexpected, due to the 
non-interacting nature of the calculation). More recently, Cebeci and Su (1988) have 
considered the full three-dimensional boundary layer problem on a prolate spheroid, 
focusing much attention on the separation aspect. Smith (1982) considered the unsteady 
counterpart of the problem considered by I, and found that a finite— time breakdown of 
the solution was possible, which he attributed to the occurence of dynamic stall. 

In an attempt to model the Cebeci at al (1980) configuration rather better, Brown 
(1985) (hereafter referred to as II) considered the problem of marginal separation of a 
boundary layer along a line of symmetry. The result of the analysis was again an integral 
equation, although this time two constants were found necessary (although these turn out 
to be related). One important assumption with this work is an a priori knowledge of the 
crossflow pressure gradient (obtained from in viscid theory), and so this work may be 
described as mildly three— dimensional, although it does have interesting and important 
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differences from the two-dimensional work of I. In particular, in general the upstream 
streamwise gradient of skin friction is no longer minus the downstream gradient (the ratio 
of these quantities being a parameter to the problem). 

The aim of this paper is to treat a fully three-dimensional marginal separation 
(again along a line of symmetry), incorporating crossflow pressure variations in the 
interaction zone. Following I and II we assume the boundary layer has a particular 
form at the separation point, which lies on the line of symmetry. Defining the 
streamwise direction to be parallel to the line of symmetry, and the crossflow direction 
perpendicular to the line of symmetry, we then make a key assumption (akin to that used 
in I and II), that the streamwise component of skin friction vanishes simultaneously (at 
the separation point) with the crossflow derivative of the crossflow component of skin 
friction. The precise details will be described in the following section. 

Work on three-dimensional interacting boundary layers (a category to which this 
work relates) has progressed somewhat slowly. Smith et al (1977) posed the 

three-dimensional triple deck problem, which they then solved for a linearised case. 
Duck and Burggraf (1986) used a pseudo— spectral method to treat the fully non-linear 
version of Smith el al (1977), and obtained solutions for a number of separated flows. 
Smith (1983) (extended by Bodonyi and Duck 1988) presented a finite— difference scheme 
for treating flows of this type, using a system of rotated coordinates, enabling the problem 
to be made quasi— two— dimensional. The amount of literature on three-dimensional 
separation remains small, in spite of its great importance; this would appear to be due to 
the great numerical difficulties encounted with flows of this type. Indeed the author is 
unaware of any detailed studies of the three-dimensional Goldstein singularity per se 
(although Cebeci et al 1980 do have some discussion of this point). 

The outcome of the analysis of this paper is a two-dimensional non-linear integral 
equation (or partial— integro— differential equation depending on the final form used). 
(The work of I and II resulted in an analogous one — dimensional equation.) In order to 
solve this equation, we shall adopt a pseudo— spectral approach, with the solution 
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technique being based partly in Fourier space, and partly in physical space. This type of 
approach is in the spirit of Duck and Burggraf (1986), although the details are very much 
different. It appears that this type of approach has considerable advantages over solution 
procedures confined to physical variables. 
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2. BASIC FORMULATION 


Let us consider an incompressible fluid of kinematic viscosity »- . 
Suppose that L denotes a typical lengthscale in the problem (for example 
the distance between a leading edge/stagnation point to the separation 
point), and that a typical freestream velocity directed along the 
centreline is U^. Then we take Cartesian coordinates (Lx.Ly.Lz), 
origin at the separation point; the flow velocity is u - U^u.v.w), 

and we define the centreline to lie along z - 0. The body surface 
(ignoring any curvature terms) lies on the plane y - 0. The pressure 
is written as pU^ p, where p is the fluid density. The Reynolds 
number (assumed large throughout this paper) is defined as 

R - IU-A. (2.1) 

The governing non-dimensional equations are then written as 

V. u - 0 (2.2) 

(u.V) u - -Vp + R-l V2u. (2.3) 

We shall suppose that all three velocity components are (generally) 
non-zero and functions of x, y, z. Because of the symmetry of the 
problem about z - 0, we expect u,v and p to be even functions about 
the centreline, whilst w wi 1 1 be an odd function of z. 

Finally, on y — 0 we require to impose the usual conditions 
of zero velocity, i.e. 

u « 0 on y-0. (2.4) 
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3. NON-INTERACTING FLOW 


As the flow approaches the point of zero skin friction (x-y-z-0) 
we suppose that the boundary layer has the standard thickness 0(R ~i). 
However, rather than allowing the boundary layer at this point to have a 
Goldstein (1948) "square-root " singularity, we assume that ambient 
conditions are such that this is avoided. This we achieve (following 1 
and 11) by a careful choice of conditions (as defined in (3.11) below), 
which results in the boundary layer exhibiting a streamwise gradient 
of streamwise wall shear (and crossflow derivative of streamwise gradient of 
crossflow wall shear) that is discontinuous at x-z-O. We now seek to 
determine the general form of this solution under these conditions, in 
this region. Taking (2.2), (2.3) and introducing classical boundary- layer 
variables , 

y - R-iy (3.1) 

v - R-iv, (3.2) 

with all other variables (both dependent and independent) remaining 
unchanged, then the three-dimensional boundary layer equations can be 
written (to leading order) 


— + 
8x 

+ 
d Y 

O 

i 

<81(8 


(3.3) 

9u 
u — 

9x 

+ V — 

av 

, du 

+ w — - 

02 

c*p o^u 

8x 8y2 

(3.4) 

dw 
u — 

0X 

+ 

3y 

9w 

+ w 

9z 

A ^0 

dp o*w 

oz 3 y^ 

(3.5) 

Sp _ 

dY 

0. 



(3.6) 

now assume the 

fol lowing 

solution expansions 

(in line with our 


earlier remarks) 
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u - U 0 (Y) + M (x2 x 2+ 2 2)J U 0 ' (Y) 

- x [wj(Y) + Vj’CY)] + \ x 2 [-V 2 ’(Y) - w 2 (Y) ] 

+ \ fi 2 (X 2 x 2 + z 2 ) U 0 "(Y) + z 2 u 23 (Y) 

- ft X (x2 x 2 + z 2)i (Vj" + wj') + ... (3.7) 

V “ v i< y ) ^ j r-— ■ Uo< y ) + xV 2 (Y) 

(X2 x 2 +Z 2)i 

-^2 X 2 x u o ' (Y) - M (X 2 x 2 +z2)i (w 1 (Y)-V 1 ’(Y)) 

X 2x2|i 

+ ; (Zwi+v^) + ... 0.8) 

(X2 x 2+ Z 2) J 

w - zwj (Y) + X z w 2 (Y) +#iz(x2 x 2+ z 2) } wi'(Y) + ... (3.9) 

P - P 0 + Pi x + i P 2 x 2 + \ P 3 z 2 + . . . (3.10) 

The key assumption here is that at X -z-0 the solution is generally 
non-di fferent iable about a small elliptic region centred about X -z-0. 

The parameter (i is essentially a lengthscale parameter, whilst 
X is a measure of the aspect ratio of the elliptic region (which may 
also be related to the relative magnitude of the streamwise and crossflow 
velocities, close to the surface). Notice that (3.7)-(3.10) satisfy 
the continuity equation (3.3). 

In order that the appropriate boundary conditions on y-0 are 
satisfied, we must have 


u o (0) - V<°> - 0 

(3.11a) 

0 

1 

0 

rH 

> 

1 

o 

'w' 

rH 

> 

(3.11b) 

w^(O) - wi'(O) - 0 (and hence Vi"(0) - 0) 

(3.11c) 

0 

1 

o 

CM 

> 

1 

O 

CM 

(3. lid) 
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V 2 * (0) - /x 2 x 2 Uo "(0) 

u 23 (0) - - 5 p 2 u o "(0)- 


(3. lie) 


(3. Ilf) 


The key assumptions here, at least from the physical point of 
view, are equations (3.11a) and (3.11c) (the remainder of these conditions 
are basically zero velocity restraints). Equation (3.11a) demands that at 
x - z - 0 the streamwise component of wall strear vanishes (c.f.I and 
II), whilst (3.11c) demands that the z derivative of cross flow wall 
shear vanishes, simul taneously, at this point. 

Since we expect the solution to be regular as Y » 0, we may 

expand the solution in the form of a Taylor series in this limit, 
for example 


®o 00 - l 

n-2 n ‘ 


*1 (Y) - l 


n-2 


n! ' 


It is a simple, matter to show that 


Y 

Vj(Y) - - U 0 f - 0^ 1 ~^L +U Q dY. 

o Uo 2 

However we expect V} to be regular as Y 
condition then also demands 


0 also, and this 


(3.12) 


(3.13) 


(3.14) 


32 “ V(°) - Pi. a 3 “ V(0) “ 0. 

*4 “ V"(0> “ 0, U 0 (0) - a 3 - 0. (3.15) 

A similar, local analysis with the z momentum equation, as 

Y r 0, reveals the solution to remain regular provided, the following 

additional conditions are fulfilled 
c 2 - w 2 ” (0) - P 3 , c 3 - 0. 

C£ - W 2 ” ’ * (0) - c 2 2 . 

(A similar analysis is also possible with V 2 and w 2 , although we 
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(3.16; 


have more than enough details to proceed to the next stage of the 
an] ays is). 

At this stage it is interesting to distinguish between the above, 

1 and II. Setting z - 0, wj - W 2 “ ... - 0, togther with /i - 1 in the 
above reduces to a system identical to that described in I. The key 
difference with the work of II is that in the latter, no restriction 
was placed on wi*(0); the primary repercussion of this was that unlike 
the above (and I), the no-slip condition downstream (x>0) was 
violated, and so an inner, Goldstein type layer was found necessary. 

This leads to a situation, in general, where the downstream wall shear is 
no longer symmetrical about x - 0 (as in the case in the above and I). 

A further difference in II was that the discontinuity in wall shear was 
centred along a line (x-0) , (rather than at a point, as above), thus 
the present structure is rather more three-dimensional. A further 
important distinction is described later. 

It is interesting to note that conditions (3.11a)-(3.11c) imply 
that at x - z - 0, all components of the vorticity vector on the wall 
are zero, and this is exactly the condition for a point of separation as 
described by Lighthill (1963). Thus the situation described above is 
important generical ly . 

Unfortunat ly , due to the limited literature available on 
integrating the three-dimensional boundary layer equations (even of a 
non-interacting nature), there do not appear to be any calculations 
illustrating the above characteristics. (Cebeci et al 1980 did present an 
example of flow past the special case of a paraboloid of zero thickness, 
with a point where U 0 * ( 0 ) - 0, but with wj'(0) + 0; this 
provided the stimulus for the work of II.) 

The discontinuity in slope of the solution here at the origin 
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(and in 1 and II) Is removed by means of a triple deck., centred at 
x-y-z-0, and we focus our attention on this aspect next. 


11 



4. THE INTERACTING REGION 


As in I and II, the appropriate streamwise scale is x-0(6), where 
6-R" 1 / 5 . We also expect the appropriate crossflow lengthscale to be the 
same, namely z-O(J), and so scaled variables (X.Z) are defined by 

X - x/6 (4.1) 

Z - z/«. (4.2) 

Consider first the main deck, comparable in thickness to the classical 
boundary layer described in the previous section. 

The expansion for the solution in this (Y-O(l)) region (guided 
partly by (3 . 7)-(3 . 10) , and partly by 1 and II) is 

u - U 0 (Y)+i{-X[w 1 (Y)+V 1 ’(Y)]+A(X,Z)U 0 '(Y)} + 5 2 {- \ X 2 IV< Y ) + w 2 (Y)] 
+ | IA(X,Z) 2 + 7(Z)]U 0 "(Y) 

+ Z 2 u 23 (Y) - XA(X,Z) [Vj"(Y) + w 1 1 (Y)]} + . . . (4.3) 

V - V X (Y) - A X U 0 (Y) 

+ b {XV 2 (Y) - AA X U 0 '(Y) + (XA) x (Vi’(Y) + w a (Y)) - (ZA) Z wj(Y)} 

+ ... (4.4) 

w - 5[Zw 1 (Y) ] + « 2 {XZ w 2 (Y) + Z A(X,Z) wj (Y)} + ... (4.5) 

P - P D + «Pj X + i b 2 P 2 X 2 + \ b 2 P 3 Z 2 + b 5 / 2 P(X,Z) + ... (4.6) 

The function y(Z) here is an arbitary function of Z, which as we 
shall see performs a key role in the investigation, and is analogous 
to the constant T in I. y(Z) is a function we may postulate, although 

we shall assume y — o(Z 2 ) as Z » a in order that a proper match is 

made with the results of the previous section. The functions A(X,Z) 
and P(X,Z) are a displacement function and pressure perturbation 
respectively, and arise from the interaction to be described shortly. 

In order that a proper match is made with ( 3 . 7 ) - (3.10) we require 
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as 


(A. 7) 


A 2 (X, 2) - /x 2 (X 2 X 2 +Z 2 ) 

The system (4.3)-(4.6) cl 
Y-0, and hence the Inclusion o 
a meaningful balance of terms 


X 2 +Z 2 -4 oo. 

early violates the no-slip condition on 
f a lower deck is necessary. In order that 
is affected, this must have thickness given by 


Y - Y6' 1 / 4 - 0(1), 


(4.8) 


(the same thickness as found in I and II). In this region, the solution 
expansion given in (4.3)-(4.5) remains essentially unchanged, except 
Y is replaced by $iY, and there is a correction to the u expansion 
of 5 2 U (X, Y, Z) , to the V expansion of « 5 / 4 V (X.Y.Z) and to the 
w expansion of 6 2 W (X,Y,Z). The governing equations for these 
correction terms are 


where 


3V 8W 

9Y 0Z 

au + 


- o, 


a 2 Y V - 

8 2 w 

av2 


dP 

0X 


a? 


a 2 u 

3y 2 


The boundary conditions to be applied to this system are that 
U, W — » 0 as Y — » ®, 


(4.9) 

(4.10) 


(4.11) 

(4.12) 


together with 

W(Y-0) - V(Y-O) - 0, 

U(Y-0) - - 22 ( A 2 + 7 <Z) - m 2 (^ 2 X 2 +Z 2 ) }. (4.13) 


Differentiating (4.10) with respect to X, and (4.11) with respect 
to Z, and adding yields 


1 * 2 y2 (£5 + ; + 

2 3X 2 dX9Z 


v dv 

>2 Y — 

ax 


a 2 p a 3 u a 3 w 

4- *f 

az 2 9 y 2 3x 8y 2 8z 


(4.14) 


a 2 p 

ax 2 
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Invoking continuity (4.9), and then differentiating with respect to Y 
yields 


3 4 V 1 

By 4 2 


a 2 Y 2 


3 3 V 

3y 2 3x 


+ a 2 


3v 

3x 


0 . 


(4.15) 


Here the problem has been reduced to a (quasi-) two-dimensional problem, 
which greatly simplifies obtaining a solution. The boundary conditions 
to be imposed are: on Y - 0 


V - 0, 

— “ 5 a 2 [2AA X - 2/* 2 X 2 X] , 

3y 


3 3 v 3 2 p 3 2 p 

3y 3 3x2 az 2 ’ 


(4.16) 


whilst V — » 0 as Y — > «. 

A very similar system was considered by Stewartson (1970), the key 
result of which shows that a solution subject to (4.16) is only possible if 


\ a 2 [A 2 + y(Z) - /£ 2 (X 2 X 2 + Z 2 )] 

£ 


<4 s > 

2^/ 2 (4! ) 


3P/3$ + J ^^2 * 

(x-Oi 


d£. 


(4.17) 


-OC 

(subject to (4.7)). 

The system is now closed by consideration of the upper deck, wherein 
Y - 0( 6 - 3 / 2 ) . The solution is standard, with the pressure being determined 


by Laplaces equation, and so following Smith et al (1977) and Duck and 
Burggraf (1986) we may write 


P (X, Z) - - ^ U 0 
The combination of 


o> 00 

(„)2 { I 3 ; A/8; 2 d US 

II K x '-S > 2 + (Z-O 2 ]* 

—00 —O0 

(4.17) and (4.18) represents a closed problem, 


( 4 . 18 ) 


the 


14 



solution of which will be discussed in the following section. 

We note that the condition (4.7), as » oc can be 

made stronger, namely 

A 2 (X,Z) - /i 2 ( X 2 X 2 +Z 2 ) -*y(Z) (4.19) 

This is analogous to the farfield condition found in I, namely 

A 2 (X) - X 2 -r (4.20) 

(where T is a constant). We note that if 7 (Z) is entirely negative, 
then we might expect there to be no difficulty in obtaining a solution to 
the system described above, with the flow remaining entirely attached. 
However as y(Z) increases, and becomes progressively more positive, (4.19) 
suggests a three-dimensional analogue to the classical Goldstein square-root 
singularity is possible, leading to a catastrophic breakdown. (Indeed, 
similar observations are noted in the work of I and II.) One of the primary 
aims of the following sections will be to determine the maximum value 
of 7 (Z) for which solutions exist. Notice, also, that as 7 (Z) 
becomes progressively more negative the effect of the interaction 
is reduced, and (4.19) becomes an asymptotic solution to the above 
system. 

In II it was assumed that the cross-stream pressure gradient 
was known a priori, (not a completly rational assumption) leading to a 
one-dimensional mathematical problem, similar to that found in I; here, 
the mathematical problem (4.17)-(4.18) is fully two-dimensional. 
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5. NUMERICAL TECHNIQUE 


At this stage it is convenient to take the double Fourier Transform of 
(4.17) and (4.18) to eliminate the pressure perturbation P (X,Z). 

De f i n i ng 

A** - J J A(X,Z)e ,kX " iCZ dXdZ (5.1) 


(with other transform quantities defined similarly), then (4.17) and (4.18) 
give 

\ a 2 [A2+ 7 (Z) -p(X 2 X 2 +z2)]** 

<-r!) 

y ( 5 )* (ik)J (k2+(>2)i U 0 («)2 a**. (5.2) 

2(|0 


This equation can be transformed back into (X,Z) space, namely 


A(X,Z) 


,7 7 L«,0 K(x) 

II [(X-0 2 +(Z-f) 2 ]i ' 


(5.3) 


where K( X ) denotes the complete elliptic integral of the first kind, 
argument x. where 


X " ^ 1 + [(X-0 2 +(Z-f) 2 ] i 1 

and 

L(X,Z)- \ a 2 [a2+ t (Z)- / x2 ( x2 X 2 + z2)], 

and 

A - 2—, 

T 3/2 (_ii) u 0 (»)2 


(5.4) 

(5.5) 

(5.6) 


However from a numerical point of view (5.3) would present a number 
of difficulties, primarily associated with singularities in the kernel 
term. Instead we prefer to adopt an alternative scheme based on the 
transform equation (5.2). If we write 
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(A,7,x,z,k,c) - ((v/sA.eV/sT.e/^/sx.e/rVsz.e-v/sk.rV/seHs.?) 


(-s'.) (f)f ,v 5 

0 - [— T- «° 2 <->] • 

( 4 !) a 2 


(5.8) 


then (5.2) simplifies to 


[ A 2 + 7(Z) - (X 2 X 2 +Z 2 )]** 
- ( i k) i (k 2 +0 2 )i A** . 


(5.9) 


In its present form (5.9) is unsuitable for numerical manipulations, 

A A A A 

because of unboundedness as X 2 +Z 2 — ;> ® (this is mirrored in k,G space 

A A 

by a singularity as k and/or G » 0). We do, however have the 


condi t ion 


A 2 - X 2 X 2 + Z 2 - 7 (Z) as X 2 + Z 2 


(5.10) 


Guided by this we write 


A - (X 2 X 2 +Z 2 +c* 0 2 )i 


7(Z)+a 0 2 
2 ( X 2 X 2 +Z 2 +a 0 2 ) i 


(5.11) 


where now we have A * 0 as X 2 +Z 2 » «. We have introduced the 

(artificial) constant a Q into the problem in order to ensure boundedness 

A A A A 

as X 2 +Z 2 » 0. (This is mirrored in k, C space by appropriate decay as 

k 2 +C 2 > oc; setting q q - 0 would cause difficulties with the transform 

in this limit.) a Q may take on any real value, and is thus a parameter 
entirely of our choice, the exact solution for A being completely 
independent of q g . 

Substituting (5.11) into (5.9) yields 


7(Z)+a p 2 + 2 ^ r + a 2 


*7> 

a[7(Z)+q q 2 ] 

! 


(ik)i (k 2 +G 2 ) i 


7(Z)+q Q 2 
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Here, for convenience we have written 

r - (X 2 X 2 + Z 2 + a c 2 )i. (5.13) 

Unfortunately (5.13) is still not quite suitable for numerical 
treatment (because of difficulties with the second term on the left-hand- 
side of (5.13) as X 2 + Z 2 > °°) . This may be rectified however by writing 


A - A 0 + B, 


(5.14) 


where 


[2A 0 r]** - ( ik) i (k 2 +C 2 )i 


f7(Z)+c*n 2 ] 


** 


(5.15) 


Details of the evaluation of A 0 are given in the Appendix, but 

A A 

now we may regard A 0 as a known function of X and Z (or A 0 ** as a 

A 

known function of k and G), and B the function to be determined. 
Equation (5.12) may now instead be written 

Jcit 

[ 7 (Z)+a g 2 ] 2 + 2B2r _ .|-y (Z )+a 0 2 j [ ^ + | ] + A c 2 + 2A 0 B + B 2 ] 

4r* J 


- (ik)i (k 2 +C 2 ) i (A 0 +B)**. (5.16) 

A A 

Equation (5.16) is cast in k,C space; however although the right- 
hand-side of this equation is most simply evaluated in its present form, in 
Fourier space, the left -hand-side is best evaluated in physical space. 
Consequently we choose to adopt a solution procedure which exploits the 
favourable aspects of both spaces. More specifically, the technique is 
as fol lows . 

A A 

For a given y(Z ) , A c , then A c ** and the inverse of 
( ik ) i (k 2 +C 2 )i A 0 ** are all evaluated (see Appendix for details). Then 

A A A A 

some "guessed" B(X,Z) (on a discrete set of X and Z points) is 
transformed into k,G space to yield B** (k.G). This is achieved using the 
Fast Fourier Transform (F.F.T.) technique of Cooley and Tukey (1965). The 
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term 


(ik)i (k2+l?2)i b** is then constructed, and transformed back to 
physical space (again by means of the F.F.T. technique). This yields an 

^ A 

estimate for the right -hand-side of (5.16) on the X, Z grid, and provides 

A A 

a means of obtaining a new estimate for B(X,Z). A number of different 
schemes were ultilised for this, but the most reliable, generally, (although 
not always the scheme that gave the most rapid convergence) was to write 

A A 

(5.16) at each X,Z location as 

B - R - qj B 2 , (5.17) 

and to use this to provide a new estimate for B, treating the right-hand- 
side of this equation as known (some under-relaxation was applied if 
necessary). In certain regions of parameter space, this iteration 

A A 

procedure was modified slightly by using (5.17) at all X, Z stations, 

A A A A 

except at one (prescribed) station, (Xq.Zq) say, B(X 0 ,Z 0 ) was 

A A 

prescribed and 7 (Z 0 ) was treated as the unknown (a similar tactic was 
followed in I). The process was repeated until the maximum change to 
the solution at any of the X, Z points fell below some predetermined 
tolerance (10“^). Since we are concerned with flows which are symmetric 
with respect to Z - 0, we need only consider Z < 0. Further, notice 

A. A. 

that if we use a gridsize in X of AX, with N x points, the grid in 
k space is given by 

Ak - 2x/NxAX. (5.18) 

This relationship arises through the use of the F.F.T. A similar 
relationship applies in the crossflow direction, 

AC - 2t/N z AZ (5.19) 

where AZ is the gridsize in the Z direction, N z the number of grid 

A A 

points, and AC the gridsize in C space. 

In the following section we go on to describe a number of results. 
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6. RESULTS AND DISCUSSION 


As a check of the general numerical technique, the two-dimensional case 
as considered in I was recomputed, using basically the same technique as 
described in the previous section. This involved solving (5.9) with Z « 0 

A A A 

C - 0, X - 1, 7 <Z) - T (a constant), with just a single Fourier 

Transform (in the X direction). Additionally, the particular details, 
given in the Appendix are rather different, although the overall scheme was 
simi lar . 

The two-dimensional case was computed with N x - 64 , AX - 0.317, and 

A A 

a D - 3; the variation of A(X-O) with T is shown in Fig.l. The value 

A 

of A here differs by a multiplicative factor t 1 / 2 from that of I, 
and the value of T here by a multiplicative factor t 2 / 5 from that of 
I (on account of the slightly different normalisation used); to facilitate 
comparison, Fig.l also indicates the values corresponding to the 
normalisation of 1 (these are denoted by a subscript I). The present 
results reproduce those of I, to within graphical accuracy. 

The non-uniqueness of the solution over a range of T, together with the 
non-existance of solution for F > 1.74 is immediately apparent. On the 
"lower branch" of the curve, computations became progressively more awkward 
as T was reduced, even using the alternative procedure where T 
itself is determined during the course of the iteration (although no 
really exhaustive attempt was made to extend these results further, the 
accuracy of the scheme having been confirmed). 

We now turn to consider a number of results for the three-dimensional 
case. To be specific, we shall focus attention on functions y{2) of the 
class 

✓s 

** A _ 7 2 

7 (Z) - T e , (6.1) 

and seel to determine the variation of soiutions with T. (Note that 
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even though this choice of *y(Z) > 0 as Z > the results will 

nonetheless be fully three-dimensional as seen by the form of the far 
field conditions (4.19)). 

/V A A A A 

Variations of A(X,Z) (at selected X,Z stations) are shown in Fig. 2. 

A A 

These results were determined with the grid N z - N x - 64, AX - AZ - 0.317, 
and a Q - 3, for the particular case X - 1. These distributions are 
qual i tat irely similar to those obtained for the two-dimensional case, 
although the range of T for which solutions exist is much increased. To 
assess the effect of numerical grid changes, Fig. 3. shows results for 

A A A A A 

A(X-0,Z-0), obtained using : a 0 - 3, N x - N z - 32, AX - AZ - 0.645 

a a 

(broken line), a c - 2, N x - N z - 64, AX - AZ - 0.317 (solidus), and a 0 - 3, 

A At 

N x - N z - 16, AX - AZ - 0.667 (dot -dashed line). 

It is heartening that there appears to be no difference 
(to within graphical accuracy) between the a c - 3 and corresponding 
a 0 - 2 results (indeed the computation was also carried out with a Q - 2, 

N x — N z * 32, AX — AZ - 0.645, and this too gave results 

indist inguishable on the graphical scale used with the corresponding q 0 - 3 
calculation). On the whole the computations agree very well, although the 
most prominant discrepancy occurs in the "nose" region of the curve, as 
may be anticipated, although elsewhere the agreement is to within the 
accuracy of the figure. The computations on the "lower branch" became 

progressively more difficult as T was reduced (a typical computation on 
the finest, most extensive grid taking of the order of ten minutes on a Cray 

XMP) , and so the curves are not extended as far as the corresponding 

two-dimensional case (Fig.l). 

Figure 4 shows the contours of constant A(X,Z) for T - 2.8031C3 
(upper branch) and for comparison Fig. 5 shows the contours for the same 

value of T, but corresponding to the lower branch (both sets of results 

~ 

were obtained using the same grid as used for Fig. 2). Since A(X,Z) 
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A A 


is proportional to the streamwise wall shear, negative A(X,Z) implies 
reversed flow. Comparing these two figures reveals the lower branch to 
involve a much larger region of reversed flow. Indeed, as the lower 
branch is followed, with T reducing, the region of reversed flow is 
expected to increase without bound (a feature observed in previous 

S' S' 

two-dimensional studies). As increases the contours became 

progressively more circular in both cases. 

A S' S' 

Figure 6 shows the variation of A(X - 0, Z - 0) with T for the 
case X - 2 (obtained using the grid N x - N 2 - 64, AX - AZ - 0.317, and 
0t o - 3). This case corresponds to a basic flow with a rather stronger 
streamwise flow near the wall compared with the crossflow. In this case 
the range of T over which the solution exists is increased over the X - 1 
results, as is the range of A(0,0). Figure 7. shows contours of constant 

A A A A A 

A(X,Z) for T - 5.342597 (lower branch) . The contours as X^+Z^ > oc i n 

this case will no longer approach the circular form of X - 1, but rather 
will tend to ellipses. 

These results presented here clearly illustrate the occurrence of fully 
three-dimensional solutions involving below - marginal skin friction. 

The basic solution technique introduced here, involving the usage of a 
pseudo-spect ral method could fairly easily be extended to other classes of 
flow, for example flows away from lines of symmetry, supersonic three- 
dimensional flows or flows where U D * (0) and Wj*(0) do not vanish 
simultaneously (such as that studied in II). However the particular 
details required in obtaining a well behaved "perturbat i on" solution 
to compute (c.f. the Appendix) would necessarily have to be modified. 
Further, and importantly, the theory presented here is not subject to the 
limitations imposed as a result of the essentially two-dimensional 
interaction condition employed in II . 
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It is worth making a few final comments regarding the choice and 

A A A A 

determination of the function 7 (Z) . In this paper 7 (Z) represents 
a quite natural extension of the constant F found in the 
two-dimensional work of I (and T+ in II). In I, T was interpreted 
as being proportional to the excess of the angle of incidence 
(or camber) over the critical value. As described in II, the 
interpretation of the two related constants r ± occuring in that study 
is less obvious. In the present work, we may interpret 7 (Z) as being 
a measure of the excess angle over the critical centreline angle; in the 
present context the local angle can depend on spanwise location, and so 
this excess will in general be a function of lateral position, which from 
a physical point of view can be regarded as a spanwise (i.e. Z) 
variation of body surface curvature and/or of streamwise conditions (we 
must however restrict spanwise variations to be on a scale z - 0(6) 
or longer). 

A A 

Our particular choice of 7 (Z) , namely (6.1) leads to necessarily 
a model problem; however because of the qualitative similarity between the 
results for this particular choice of 7 (Z) and those found in I and 
II, we fully expect the present three-dimensional results to give a good 
indication of the qualitative feature of general flows of this class. 

A A 

Since 7 (Z) represents an "excess" over the critical state, 

by its nature it controls the amount of interaction in the problem (see 

also the remarks, towards the end of Section 4 regarding this point). 

Hence it is not directly deducible from "non-interacting" calculations such 
as those of Cebeci at al (1988). (The constant T found in I is also 
not directly obtainable from "non-interacting" calculations.) 


23 



Acknowledgement s 

The author wishes to thank Professor O.R. Burggraf for encouragement 
and many useful discussions relating to this work. Computations were 
carried out at the Ohio Supercomputer Center, the University of Manchester 
Regional Computer Centre (with computer time provided under S.E.R.C Grant 
No. GR/E/25702) and at NASA Lewis Research Center (whilst the author was 
visiting under the I COMP program). 


I 


24 


Appendix 


Here we consider the evaluation of A 0 and A 0 ** . Since certain of 
the components in the evaluation of A 0 are unbounded as x2+z2 » re 

A A 

(and this is reflected by singularities as k2+£>2 > 0 in Fourier space), 

a certain amount of analytic work is necessary, prior to any numerics. 
Fortunately, the inclusion of the a Q term (see (5.11)) ensures 

decay in k, C space as k^+C^ > and so this limit requires no 

special at tent ion. 

By (5.15) we have 


(rA 0 )** - i (ik)i (k2+fi2)i j r _ _I T 1 . j 


(A . 1 ) 


It is straightforward to show 

r ** - _ 2ir (k2 + c2 x 2 ) -3/2 x 2 [j + ^-^2+^2)* 

(assuming a D > 0) . (A. 2) 

The quantity r** (ik)i (k 2 +fi2)J ls transformed back to X,Z space 

by the integral 

2 °? x/2 . vi a n p2cj) 

J J v c °ssp; — e x (i + a 0 p 24 ,) cos (p2 C os^X+£) cos (p2sin^Z)dpdp 


o o 


* 3 X 

(A. 3) 

where $ - [cos2^pf\2sin 2 y>]*. (A. A) 

A A 

Here the Fourier variables k, C have been transformed to Fourier polar 
variables p , <p t and then the radial variable p is introduced by 
p * p* . The result, (A. 3) is a fairly well behaved integral, which is 
readily amenable to numerical quadrature (the weak square root 
singularity at ~ x/2 was not found to cause any difficulties). 

The n a 0 2/2r M term in (A.l) was treated in much the same way as 
above, using the result 


( 


*,« - 2,(iVA2)-i e -“° 


(A. 5) 
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The remaining term was evaluated rather more numerically. We have 



[ ]** “ 2 I 7 (Z > cosCZ K 0 

o 

1 I <z 2 ««„ 2 >i] 

| dZ 


(A. 6) 

(having 

carried out the X transform. 

Here K 0 ( X ) 

denotes 

the 

mod i f i ed 

Bessel 

function of order zero, argument The Z 

int egrat ion 

genera 1 ly 

requires a numerical treatment (depending on the particular 

form 

for 

7<Z)). 

Once evaluated, the quantity 

(ik)i(k 2 +0 2 )* 

f 7< z > ] 

l 2r J 

was then 


transformed to (X,Z) numerically, completing the determination of A 0 . 
A 0 ** was then evaluated without any special numerical treatment, 

A A A 

( ik) i (k 2 +£? 2 ) } A q ** was constructed, and then transformed back to X,Z 
space (again, with no special treatment), thus completing the preliminary 
computation, prior to the iteration commencing. 
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FIGURE 1. - VARIATION OF A(0) WITH T (TWO-DIMENSIONAL CASE). A,, f, DENOTE VALUES 
USING NORMALIZATION OF I. 
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FIGURE 4. - CONTOURS FOR CONSTANT A(X,Z) y(Z) = 2.803103 e' Z , A= 1, UPPER BRANCH SOLUTION. 



Z 

- - - - - _ 7 2 

FIGURE 5. - CONTOURS FOR CONSTANT A(X.Z), y(Z> = 2.803103 e 1 , LOWER BRANCH SOLUTION. 
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FIGURE 7. - CONTOURS FOR fl(X = 0, Z = 0). y(Z) = 5.342597 e* Z A = 2, LOWER BRANCH SOLUTION. 
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